Draft version March 2, 2013 

Preprint typeset using I^TgX style emulateapj v. 10/09/06 



O 

< 

oo 



Oh 

6 



(N 
> 

O 

in 

m 

^' 

o 

(N 



X 



MAPPING EARTH ANALOGS FROM PHOTOMETRIC VARIABILITY: 
SPIN-ORBIT TOMOGRAPHY FOR PLANETS IN INCLINED ORBITS 

YUKA FUJII^ AND Hajime Kawahara^ 
Draft version March 2, 2013 

ABSTRACT 

Aim ing at obtaining detailed information of surface environment of Earth analogs jKawahara fc Fujiil 
(|2Qllf ) proposed an inversion technique of annual scattered light curves named the spin-orbit tomog- 
raphy (SOT), which enables us to sketch a two-dimensional albedo map from annual variation of the 
disk-integrated scattered light, and demonstrated the method with a planet in a face-on orbit. We 
extend it to be applicable to general geometric configurations, including low-obliquity planets like the 
Earth in inclined orbits. We simulate light curves of the Earth in an inclined orbit in three photomet- 
ric bands (0.4— 0.5/im, 0.6— 0.7/im, and 0.8— 0.9/im) and show that the distribution of clouds, snow, 
and continents is retrieved with the aid of the SOT. We also demonstrate the SOT by applying to 
an upright Earth, a tidally locked Earth, and Earth analogs with ancient continental configurations. 
The inversion is model-independent in the sense that we do not assume specific albedo models when 
mapping the surface, and hence applicable in principle to any kind of inhomogeneity. This method 
can potentially serve as a unique tool to investigate the exohabitats/exoclimes of Earth analogs. 

Subject headings: astrobiology - Earth - scattering - techniques: photometric 



1. INTRODUCTION 

One of the ultimate goals in exoplanet study is the dis- 
covery of extraterrestrial life. Recently, the abundance 
of Earth-sized exoplanets has be en implied from th e 
radial velocity observations f e.ps;.. iHoward etaTI 2010l ). 
the transit observations (e.g.. iBorucki et aD'2 011), and 
gravitational planetary micro- lensing (e.g., Cas san et al.l 
[2Q12t) and some in so-called habitable zones (HZs) have 
been already detected. In near future, direct detection 
of planetary light will play a crucial role in exploring 
the surface environment and meteorology of planets in 
habitable zones (HZ planets), and in detecting possi- 
ble biosignatures there in cluding O2, O3, H2 O absorp- 
tion lines (e.g., Owen 1 98 0: 1 Angel et aL|| )1986: Leger et al.l 
I1993I : ISelsis et al.l 120021 : iKaltenegger et al.l l2010|, and 
references th erein) and the vegetati on's red e dgefe.g., 
Arnold et al.ll2002l: IWoolf et al. 2002; Seager et al. 2005; 
Montahes- Rodriguez et al. [2OO6I : Hamdani et al.i i2006i : 
Kiang et al. 2007a.b). 

As a test bed for future observations of HZ planets, 
[Ford et al.. (2001) simulated the diurnal variation of the 
scattered light of the Earth. They found 20% (in re- 
alistically cloudy cases) or 150% (in cloud- free cases) 
variations originated from the inhomogeneity of surface 
components and cloud cover. Their results encourage 
the possibility to reconstruct the landsc ape of exoplan- 
ets fro m scattered light curves. Later on. lOaklev fc CasE 
(|2009 h suggested the inversion using the ocean glint at 
crescent phase and the albedo difference between ocean 
and land. Their methods work reasonably well when the 
cloud cover is lower than the actual Earth. ICowan et al.l 
(|2009l . l201lf ) performed a principal component analysis 
(PC A) on the multi-band photometry of the Earth ob- 
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tained with EPOXI mission. They found that a domi- 
nant eigencolor of the Earth corresponds to the cloud- 
free land component, which allows one to draw a one- 
dimensiona l map of c ontinental distribution along the 
longitude. Fujii et al.l (|2010[ ^2011) decomposed of the to- 
tal scattered light of the Earth into several albedo models 
— ocean, snow, soil, vegetation, and clouds — and found 
that the extracted time variation of them roughly re- 
covers the actual surface inhomogeneities including the 
distribution of vegetation as well as continents/ocean. 

Recently, Kawahara & Fujii (2011, hereafter Paper I) 
proposed a new inversion method named "spin— orbit to- 
mography" or SOT, to sketch a two-dimensional map 
of the exoplanet 's surface using both spin rotation and 
orbital revolution of the planet. This method directly 
maps the reflectivity over the surface with the Tikhonov 
regularization, which balances between the observational 
noise and spatial r esolution of the su r face. This was an 
improvement over iKawahara fc Fujiil (|2010l ) , where the 
inversion was regularized by the bounded variables least 
square method and thereby speciflc albedo models were 
required. We stress that the SOT retrieves the spatially- 
resolved image of exoplanets from the disk-integrated 
light curve without directly resolving of the planet. In 
this sense, the SOT can be a precursor to the projects 
which attempt to resolve planetary image directly with 
very long baseline, such as Hyperteles cope (150km ar - 
ray size with 150 apertures of 3 m; e.g. lLabeyri^ll999f ). 
In a hypothetical case where a mock Earth with plane- 
tary obliquity 90° is in a face-on orbit, they found that 
the albedo map reconstructed from a single photometric 
band traces the distribution of clouds. They also showed 
that cloud-free surface features could be recovered from 
the difference of two photometric bands since the cloud 
reflection spectra are roug hly constant over optical wave- 
lengths (e.g. Figure 8 of lFuiii et ani201lh . The plane- 
tary obliquity was also reasonably estimated within the 
framework of the SOT. 

This paper aims to further develop the SOT by con- 
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sidering various configurations, including low-obliquity 
planets like the Earth in inclined orbits. In the case 
where the planetary orbit deviates from face-on, the 
phase angle of the planet is changed and the anisotropy 
of scattering, primarily by clouds and the atmosphere, 
can have non- negligible effects. In this paper, we take 
account of such effects to make the methodology appli- 
cable to more general cases. 

The organization of this paper is as follows. In Section 
[21 we define the geometric configuration and describe the 
methodology of the SOT in detail. In Section [3l we pre- 
pare mock light curves to be applied to the demonstra- 
tion of the SOT. The simulation scheme, adopted data 
sets, and examples of simulated light curves are found 
there. Section |4] shows the results of two-dimensional 
mapping and obliquity measurements via the SOT, and 
estimates instrumental requirements for future missions. 
Case studies with different configurations are demon- 
strated in Section [5] with examples of a zero-obliquity 
planet in more inclined orbit and a tidally-locked planet 
(Section l5.ip as well as Earth models with ancient conti- 
nental distribution (Section [52]) • Finally, Section [6] sum- 
marizes the results of this paper and discusses the impli- 
cations of the SOT. 

2. SPIN-ORBIT TOMOGRAPHY 

In this section we describe the methodology of the 
SOT, which was originally proposed in Paper I. While 
we described essentials of the SOT in Paper I assuming 
the face-on orbit, in this paper we extend our method 
to be applicable for an arbitrary direction of the line of 
sight. This generalization adds the orbital inclination i 
to the set of geometrical parameters. 

2.1. Geometry 

Figure [1] illustrates the geometrical configurations. We 
define the spherical coordinate system, 6 (colatitude) and 
(j) (longitude), fixed on the planetary surface as shown 
in Figure [T] (a). The spin motion is specified by <I>(t), 
which is defined as the angle between the direction of 
the equinox and the origin of (j). Denoting the angular 
velocity of the spin rotation from observation by cjgpin, 
we can write ^{t) = cjspin^+ ^offset- Note that <l>offset only 
changes the prime meridian and have no physical conse- 
quence. We will discuss the measurement of the spin ro- 
tation by the autocorrelation in Section 14. li The orbital 
phase is described by O, taking the superior conjunc- 
tion as its origin. We adopt two parameters to specify 
the planetary spin axis: the obliquity, and the orbital 
phase of the vernal equinox, Oeq. The inclination of the 
orbit is indicated by i. 

Throughout this paper, we assume that the planetary 
orbit is already determined, i.e., inclination z, the time 
revolution of O, and the phase of superior conjunction are 
known. Geometrical parameterizations are summarized 
in Table [H 

We introduce three unit vectors to describe the direc- 
tion of the incident ray and the observer: the vector from 
the planetary surface to the host star, es, the vector from 
the surface to the observer, eo, and the vector from 
a planetary center to the surface, cr as shown in Fig- 
ure [T] (a). For mathematical simplicity, we can take an 
auxiliary coordinates where the x-axis is the direction of 
equinox and the z-axis is the orbital axis (Figure [T] (b)). 



TABLE 1 

Coordinate System and Geometrical Configuration 

Symbol Explanation 
(Coordinates) 

Colatitude fixed on the planetary surface 
(f) Longitude fixed on the planetary surface 

B Orbital phase measured from superior conjunction 
<l> Phase of spin rotation 

(Configuration of the System) 

Planetary obliquity 

1 Orbital inclination 

Oeq Orbital phase of vernal equinox 

and the components of the three unit vectors are: 

es = (cos(e - Oeq), sin(e - O,^), Of, (1) 
eo = (sin i cos Oeq, — sin i sin 0eq, cos i)'^ , (2) 




cos ((/) + $) sin 6> \ 
cos^sin (0 + sin6> + sinCcos^ , (3) 
— sin ( sin {(j) + sin + cos ( cos J 



where R{Q is the clockwise rotation matrix around x- 
axis and e'^{(j)^0) = (cos sin ^, sine/) sin ^, cos^)^. 

2.2. Modeling of Scattered Light 

Scattering property of each patch of the planetary sur- 
face is generally expressed by the bidirectional reflectance 
distribution function (BRDF), /('^^o, i^i, which is a 
function of the incident zenith angle i^o, the scattering 
zenith angle i^i, and the relative azimuthal angle between 
the incident and scattered light cp (if the surface have no 
special direction) as shown in Figure [T] (a). The scattered 
intensity at each planetary pixel is expressed as 

^ =i^*£i^'/K^o,^i,(/^)cosi?ocos^i, (4) 

where /^(i^o, '^^i, ^) is the BRDF at wavelength F^£ 
is the incident flux at Rp is a planetary radius, and 
duj = sin 6d0d(j). Note that i^o, and (p are functions 
of all of the geometrical parameters listed in Table [TJ 

For the reflection model, we assume that the scatter- 
ing is isotropic (Lambertian) . Then the BRDF becomes 
independent on i9o , 'i^i, or and 

/^(LAM)(^0,^1,^) = ^^^. (5) 

where ai{(j)^0) is the albedo at the planetary surface at 
wavelength I. Under the above assumption. Equation (|4]) 
reduces to 

Ii = ^^ J a,{^,e)WiWvduj, (6) 

where the weight functions for the illuminated and visible 
area are defined by Wj = maxjco sT^n, Q| g^nd Wy = 
maxjcos'i^i, 1}, respectively (see also lCowan et al."2QQ9^. 

Using the three unit vectors (Equations ([I]) — ®), we 
obtain 

Wi{(^, 0; ^, 6) = max{es • br, 0}, (7) 
V^y (^, 6>; ^) = max{eo • br, 0}. (8) 

The explicit forms of Equations ([7]) and (|8|) are given in 
Appendix A. 




Fig. 1. — Schematic configurations of the planetary system and surface. Panel (a) shows the spherical coordinate fixed on the planetary 
surface: the complementary latitude, 9^ and the longitude, (j). The spin rotation axis is indicated by ^ = 0. We denote the spin motion by 
<l>, the angle between (0 = 0,^ = 7r/2) and the direction of equinox. We define three unit vectors on the planetary surface (the normal 
vector, Ci^, the vector to the observer, cq, and that to the host star, es) which specify the arguments of the BRDF /(i?o, 't^i, V^; the 
incident zenith angle, i?o, the zenith angle of scattering, and the relative azimuthal angle between the incident and scattered light, 
Panel (b) illustrates the geometrical configuration in the orbital frame. The orbit lies in the x — y plane, corresponding to the directions of 
the equinox (x-axis) and the solstice (y-axis) of the planet. The 2;-axis is defined by the normal of the orbital plane. The obliquity C, and 
the inclination i are defined as the angle from the 2;-axis to the spin axis and to the line of sight, respectively. The orbital phase measured 
from the superior conjunction (SC) is denoted by B. The orbital phase of equinox is indicated by Beq- Panel (c) displays the geometrical 
configuration in the observer frame where the longitude of the ascending node VL is indicated. 

We assume that the pixel size is smah enough to satisfy 



We define the weight function for the ihuminated and 
the visible area as follows: 

W{^, 0; 6) = Wi{^, 0; 9) lyy ((/), 0; (9) 

With this weight function, Equation (j6j) reduces to 

h{^, 6) = / W{^, 6; ^, 6)a^((/), e)duj. (10) 

The integral of the weight function over the planetary 
surface provides the phase function of a uniform Lambert 
sphere (e.g., Russell 1906), 



LAM 



(„,=/ 



= / Wduj 



= - [sin a + (tt — a) cos a] , (11) 
o 

where a is the phase angle, defined as 

cos a = es - eo = s'mi cos 6. (12) 

Figure [2] displays the behavior of the weight func- 
tion in the case of (C = 90°,eeq = 0°,i = 0°), (C = 
23°.45,eeq = 90°, z = 45°), and (0°, 90°, 60°). The 
weight function covers the region of < 6 < Oy,- (or 
6>y,+ < 6> < TT ) and -180° < < 180° as the planet ro- 
tates around the host star and its spin axis, where Ov,± 
is the boundary of Wy described in Appendix A. 

2.3. Solving Inverse Problem 

In order to linearize Equation (p!Q|) . we discretize the 
surface to Mpix small pixels. The center and the area of 
the jth pixel are denoted by {(j)j^Oj) and 5j, respectively. 
Equation ([10]) is then expressed as 



W{(l),e;e;^)duj. 



W{<l>,e;e,<^)du 



« F^eRl V ^^hMw{^„0^;e, *)A^„ (14) 

^ — ^ TT 

j 

where Auj denotes the solid angle of jth pixel. 

We consider A^data bins of observational data. We de- 
fine the i-th data point {i = 1,2,..., A'data) as 



(15) 



where Si is noise of the i-th data point. 

After these discretizations, we obtain a hnear discrete 
equation 



d = Gm + e, 

Gij = W{(j)j,ej]Qi,^i)l^ojj, 
mj=Rlae{<pj,0j). 



(16) 
(17) 
(18) 



(13) 



As a pixelization scheme for the planetary surface. Hi- 
erarchical Equal Area isoLatitude Pixelization (HealPix) 
is adopted in this paper^. Unless otherwise noted, the 
resolution parameter for HealPix A^side is set to 8 corre- 
sponding pixel numbers of which is Mpix = '^'^^side ~ 
768 (Gorski et al. 2005) . 

Note that when the orbit is not circular one needs to 
caution about the incident flux F^£ because it changes 
in time. However, the procedure of analysis after trans- 
lating the observed planetary intensity into d and the 
observing time into {O, ^} remains essentially same. 

2.3.1. Tikhonov Regularization 
^ http:/ /healpix.jpl. nasa.gov/index.shtml 
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a) ^ = 90°, 0eq= 0°, / = 0° 



b) ^ = 23.45°, 0eq= 90°, / = 45° 



C) ^ = 0°, 0eq= 90°, / = 60° 
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Fig. 2. — Examples of geometrical dependence of the weight function VK(0, 0\ ^\ B). Annual variations are indicated by different thin lines. 
Solid, large-dashed, middle-dashed and small-dashed lines indicate B = 0°, 90°, 180°, and 270°, respectively. Inner and outer contours 
indicate VK((/), ^; <l>; B) = 0.05,0.3 and 0.6. Thick lines indicate the boundary of the visible area, d = Ov,±{4>) in Equation (|A4|) . Panel (a) 
indicates an extreme case: (C = 90°,Beq = 0°,i = 0°), corresponding to the geometry considered in Paper I. Panel (b) corresponds to 



the fiducial case of the Earth-twin we consider in this paper: = 23°. 45, Be 



on (/) + <^», in other words, longitude minus the phase of the spin rotation, we adopt 

showing diurnal variations explicitly. Panel (c) indicates the upright case or the tidally locked case: (C ■ 
tidally locked case, the a;— axis becomes and does not depend on 

For the time being, we assume that all geometrical pa- 
rameters (i.e. not only i but also ( and 6eq) are known. 
To solve the linear inverse problem (Equation (p!6|) ). we 
add a regularization term to the misfit function; the reg- 
ularized solution is given by minimizing 



: 270°, i = 45°). Since the weight function solely depends 
to the x— axes in panels (a) and (b) instead of 
: 0°,Beq = 270°,i = 60°). In the 



data on estimation o f rriest.A^ we int roduce the integrated 
sensitivity vector S (|Zhdanovll2QQ2f ) defined as 



ATdata 

Q= E 

i=l 



Mi - Gi 



(19) 



where cr| is variance of the observational error and P{m) 
is the penalty function, which measures the regularity of 
the solution. The regularization parameter A controls 
effective degree of freedom of the model. 

In this paper, we adopt the Tikhonov regularization, 
also known as the du mped least-square method (e.g., 
lMenkelll989l : lTarantolal[2QQ5l : lHansenll2QlQf ). The penalty 
function of the Tikhonov regularization is: 



^tik(m) = |m - m|' 



(20) 



where rfi is the prior of the model, each component of 
which is fixe d to the average of apparent albedo (e.g., 
iQiu e t alT 2003) times R^. The solution which minimizes 
Equation (fT9|) with Equation (|2Q]) is given by 



^est,A 



■ Grfi) + rh. 



(21) 
(22) 



where we define di = di/ai and Gij = Gij/cii. The or- 
thogonal matrices U and V are given by the singular 
value decomposition of G = UAV^ and Ki is the ith 
eigenvalue of G, tha t is, the i-th. component of the diag- 
onal matrix A fe.g.. ' Menkelll989l : lHansenll2010f ). The 6ij 
is the Kronecker delta. 

If data have enough information to reconstruct the 
model, rfi has little infiuence on mest,A7 while if data 
have almost no information, mest,A converges with rfi. 
From the Bayesian viewpoint, the Tikhonov regulariza- 
tion is equivalent to the parameter estimation with a 
Gaussian-type prior with the average rfi and the covari- 
ance matrix We briefly summarize the statistics 

of the Tikhonov regularization in Appendix C. 

The probability that each pixel illuminates depends 
on the obliquity, the direction of the observer, and the 
latitude of the pixel. During a year of the planetary 
system, some pixels are frequently sampled and others 
are not. In order to quantify the relative influence of 



(23) 



The estimated value of rrij is sensitive to data if Sj is 
large and vice versa. Clearly, mest,A,j should be equal to 
771 j if Sj is 0. 

The regularization parameter A is chosen based on 
"the L-curve criterion", which is the maximum curva- 
ture point of the model norm |mest,A —ffil versus residu- 
als \d— Gmest,x\ plot (|Hansenll2010[ ). The details of the 
L-curve criterion is described in Appendix D. 

2.4. Obliquity and Equinox 

So far we have described the methodology to recon- 
struct the albedo maps with known geometric param- 
eters. In reality, the weight function W and thus the 
design matrix G are functions of ( and 6eq as well: 

d = G{C,Gec)m^s, (24) 
G,,(C, Oeq) = W{^j,Oj',G,, C, 6eq)Ac^,. (25) 

While direct imaging will be able to estimate i and the 
superior conjunction, two parameters which specifles spin 
axis (obliquity, and orbital phase of equinox, Beq) are 
diflicult to estimate by any other known techniques. In 
face-on cases. Paper I showed that these two parameters 
can be simultaneously estimated via the SOT itself by 
taking them as free parameters and minimizing Q^. In 
the same manner, we search for the best-flt ( and Oeq by 
the Nelder-Mead method with varying A. At each trial 
value of C and 6eq, the best-flt model vectors rrij are 
derived by Equation (|2T]) . The set of the best-flt ©eq 
and rrij for different A allows us to draw the L-curves, 
and hence we can flnally determine the appropriate A by 
the L-curve criterion. 

3. SIMULATION OF SCATTERED LIGHT CURVES OF AN 
EARTH-TWIN 

In this section, we describe our Earth model to be input 
to the SOT as well as the resultant light curves. 

^ In the case of face-on orbit discussed in Paper I, B was mea- 
sured from the vernal equinox instead of the point of superior con- 
junction because the latter cannot be defined, and accordingly the 
measurement of B^, the phase of the start of observation relative 
to the equinox, was to be estimated rather than Beq. 
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So far, various authors have carried o ut simulations of 
scattered lig:h t curves o f an Earth- twin (iFord et ari|2QQl | : 
Tinetti et all l2006a''b'; 'Montahes-Rodrfeuez et al." '2006'; 
Palle et all [2 008: Oaklev & Cash 2009; Arnold et al. 
20091 : iDou ghtv Wolfil2010M Fuiii et al. 2010, 2011). In 
thi s paper, w e follow the simulation scheme described 
bv iFujii et al.l (|201l[ ): we (1) pixelize the surface of the 
planet into 2°.0 x 2°.0, (2) assign the parameters of the 
scattering properties to each patch according to the data 
obtained with the MODIS (MODerate resolution Imag- 
ing Spectroradiometer) onboard the Earth Observing 
Satellites Terra/ Aqua, (3) ca lculate the radiative trans- 
fer in each patc h by rstarGh (|Nakajima fc Tanakal 119881 : 
iNakajinial Il983f ). and (4) integrate the emergent light 
from each patch over the illuminated and visible portion. 

Since we do not have a specific instrumental design 
in mind, we simply assume three 0.1/im width bands: 
0.4-0. 5/im (Blue), 0.6-0. 7/im (Orange), and 0.8-0.9/im 
(NIR). Since it is fairly time consuming to sum up the ra- 
diance over planetary surface and to integrate fine spec- 
trum over photometric bands, we create a lo ok-up table 
with r elevant parameters in the same way as iFujii et aU 
(|201lf ) and estimate the radiance at a given pixel by lin- 
early interpolating it. 

In Section 13. If we describe the adopted data sets and 
our assumptions in more detail. Examples of simulated 
light curves are shown in Section [321 

3.1. Assumptions and Input Data 

The rstarGh calculates the radiative transfer through 
the atmospheric layers given the properties of atmo- 
sphere, ground albedo, the direction of the host star, 
and that of the observer. The vertical profiles of pres- 
sure, temperature, and molecule composition in the at- 
mosphere are assumed to be uniform for simplicity; we 
adopt the US standard model. ^ The thermal profile af- 
fects little on the resultant spectra at the wavelengths 
where the thermal radiation is not dominant. 

Since cloud cover significantly changes the refiection 
spectra, we insert a cloud layer according to the actual 
cloud data of each pixel. In our model, cloud cover is 
specified by two parameters: cloud cover fraction /dd 
and cloud optical thickness Tdd- The daily global map 
for these two parameters is taken from the Terra/MODIS 
Atmosphere Level 3 Product^. Specifically, we adopt the 
data of 2008. Other parameters for cl oud optical proper - 
ties are fixed to the value adopted by iFuiii et al] (|201l[ ) . 
In particular, we assume that clouds are composed of 
pure liquid water and that the cloud layer is at 3.5—6.5 
km. The altitude of the clo ud layer affects th e continuum 
level little (e.g. Figure 8 of lFuiii et aD l201lV 

The boundary condition of the radiative transfer is 
determined by the BRDF of the surface. The MODIS 
project adopts the Ro ss-Li model to sp ecify the BRDF 
of land surface (e.g., ILucht et al.l |200Q[ ) and offers the 
three coefficients in the Ross-Li model (coefficients of 
isotropic, geometric, and volume terms). We, however, 
use the isotropic term only in our simulation to match 
the specification of rstarGh. We obtain the value of the 
coefficient of the isotropic term at each pixel from "snow- 

^ U.S. Standard Atmosphere, 1976, U.S. Government Printing 
Office, Washington, DC, 1976 

^ http:/ /ladsweb. nascom.nasa.gov/ 



free gap-filled MODIS BRDF Model Parameters" . These 
data sets are a spatially and temporally averaged prod- 
uct derived from the 0.05 resolution BRDF/albedo data 
(v004 MCD43C1)^. The original data sets are merged 
into the adopted 2°.0 x 2°.0 resolution by averaging. 

In order to incorporate the effect of seasonal variation 
of snow cover, we adopt the monthly snow cover map 
offered by MODIS project.^ We replace surface albedo 
by the reffection spectrum of fine snow for pixels with 
snow cover fraction larger than 50%. 

In short, our model in this paper incorporates the effect 
of time variation of the cloud cover fraction, cloud optical 
thickness, surface BRDF, and snow cover. 

3.2. Light Curves 




-90 m \m 270 

Orbital Phase © [deg] 



Fig. 3. — Yearly light curves of an Earth- twin in the case of 
(C = 23°. 45, Beq = 90°, i = 45°). Signal-to-noise ratio (S/N) 
is imposed so that S/N after folding the light curves for 6 
days becomes 20). Top: light curve of Blue (0.4— 0. 5/im) band and 
phase curves of a Lambert sphere with albedo 0.33. Middle: light 
curve of Orange (0.6— 0. 7/im) band and phase curves of a Lambert 
sphere with albedo 0.25. Bottom: light curve of NIR (0.8— 0.9/im) 
band and phase curves of a Lambert sphere with albedo 0.28. The 
inserted panels enlarge the diurnal variations at B ~ 15° and show 
theoretical light curves by solid lines for reference. In the analysis 
in Section m data within the region bracketed by dot-dashed arrows 
are used. 

Figure [3] displays the resultant light curves of the Earth 
in three photometric bands with signal-to-noise ratio 
(S/N) = 20/V6 in the case of (C = 23°.45, Geq = 90°, i = 
45°). The inserted panels show the diurnal variations and 
theoretical (noiseless) light curves in black dotted lines. 
The variation is due to the spin rotation and reflects 
the longitudinal inhomogeneity. The amplitude of varia- 
tion 20%) is con sistent with former works (|Ford et al.l 
[200TI : lOaklev L Cash,2009. ) . 

^ Available at |http://modis.gsfc.nasa.gQv/| 
^ |http: / /modis-snow-ice. gsfc.nasa.gov/ 1 
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Fig. 4. — Top: yearly variations of the difference between the 
reflectivity in NIR (0.8-0.9/im) and that in Blue (0.4-0.5/im), 
both of which are shown in Figure [3] The geometric parameter 
set is (C = 23°. 45, Beq = 90°, i = 45°). Bottom: same as the top 
panel but the difference between the reflectivity in NIR (0.8-0.9/im) 
and that in Orange (0.6-0.7/im). 

The yearly variation is due to the change in the ihu- 
minated and visible area (Section [2]), i.e., waxing and 
waning of the planet. In Figure [3l yearly variations of a 
hypothetical Lambert sphere are also plotted by dashed 
lines for reference. Simulated light curves substantially 
exceed the Lambert phase curves at phase O = 180°, 
corresponding to the maximum phase angle. The excess 
at this phase is likely due to the enhanc ed forward scat- 
terin g: by clouds and atmosphere (e.g., iRobinson et al.l 
|2Q1Q[ ). This deviation from isotropic scattering possibly 
causes a systematic bias on the reconstructed map. Un- 
less the change of phase angle is negligible, one must 
consider this anisotropic effect. In this paper, we avoid 
this problem in a simple way: we do not use the data at 
crescent phase since the fractional deviation from Lam- 
bert sphere is quite significant at this phase (see Section 
14.21 below for further discussion). 

Figure m shows the differential light curves between two 
photometric bands (NIR— Blue, NIR— Orange) for later 
use. 

4. GLOBAL MAPPING OF AN EARTH-TWIN 

In this section, we apply the SOT to the simulated light 
curves. In Section 14. 1[ we begin with pre-processing of 
data to be relevant to the SOT by measuring of the spin 
rotation period. A two-dimensional mapping is discussed 
in Section 14.21 The obliquity measurement is studied 
in Section l43l The desired aperture size for reasonable 
mapping via the SOT is estimated in Section l474l 

4.1. Pre-processing 

First of all, we consider translation of the observation 
time to phase angles {0,<l>}. Orbital phase O is known 
once the planetary orbit is determined, which is one of 
the assumptions in this paper. On the other hand, spin 
phase ^ needs to be known by measuring spin period 
'^spin- Therefore, we perform aut o- correlation anal ysis 
on our mock data as suggested bv iPalle et al.l (|2008l ). 

Figure [5] displays the auto-correlation coefficients from 
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Fig. 5. — Autocorrelation coefficients computed from the ffist 14 
day (red solid), 30 day (green long-dashed), and 60 day (blue short- 
dashed) mock observations in NIR in the case of ((" = 23°. 4, Beq = 
90°, i = 45°) (Figure [g. 

0.8 to 0.9 /im data obtained in the first 7 day, 30 day, 60 
day mock observations. We find that the spin rotation 
period of the Earth, 24 hr, is safely measured from any of 
single-band observations for ~ 1 month, which enable us 
to assign ^. Note that the autocorrelation coefficient at 
t 7^ 24 (hr) is not approaching as the number of days 
increases because of the phase variation. When z = 0, 
the phase angle does not change and the autocorrelation 
is conversing to 0. 

Furthermore, we stack the light curves for every 6 days 
to increase the signal-to- noise ratio per data point. After 
all we are left with 30 x (366/6) = 1830 data sets of 
{e,^,/i,cri} for each band (Ar=1830) with S/N=20. 

4.2. Two-dimensional Mapping 

We now adopt the SOT to mock light curves of the 
Earth. First, we naively used one-year light curves with- 
out considering the anisotropic effect. The resultant map 
displayed in panel (a-1) in Figure [6] exhibits a distinct in- 
tense region at the North pole and hardly traces any real 
features of the Earth. This is due to the fact that the 
reffectivity at large phase angle (O ~ 180° see also Equa- 
tion (p!2] )) is much larger than expected from Lamber- 
tian phase curves (see Figure [3]). In reality, the increase 
in reffectivity at ~ 180° originates from the forward 
scattering of clouds and atmosphere as mentioned above. 
This result implies that we cannot ignore the bias from 
the anisotropy of scattering in the case of inclined orbit. 
In contrast, the anisotropy mostly leads to the offset for 
the face-on case {i = 0) discussed in Paper I since the 
phase angle is constant. 

In order to reduce this bias, we exclude the data with 
large phase angle a. In this paper we use the data that 
satisfies the condition a < 90° where the deviation from 
the phase curve of a Lambert sphere is significant (see 
also Figure 2 of iRobinson et aITl2010f ). This condition 
corresponds to —90° < Q < 90° (the bracketed region in 
Figure [3|) based on Equation (p!2|) . 

Panel (a-2) in Figure [6] shows the albedo map recon- 
structed from the Blue light curve at a < 90°. The recon- 
structed map exhibits the highly reflecting zone at high 
latitude in the north hemisphere. The high-albedo re- 
gion at 6> ~ 60° can be interpreted as a pattern of clouds 
and snow, which have high reflectivity in the Blue band. 
We can compare the reconstructed map with an annual 
mean distribution of cloud optical thickness (panel (b)) 
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(a- 1 ) Blue (SN=20), no cut Sensitivity (b) Annual Mean of Cloud Optical Thickness 
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Fig. 6. — (a-1) Two-dimensional mapping from one-year observation in Blue band with S/N=20 in the case of ((" = 23°. 4, Beq = 90°, i = 
45°) without consideration of bias caused by anisotropic scattering, i.e., data at all phases are used in the analysis, (a-2) Same as the top 
panel but data at 90° < B < 270° are not used in the analysis since the effect of forward scattering of clouds is not negligible, (b) Annual 
mean of cloud optical thickness, (c) Cloud cover fraction in 2008 March, (d) Two-dimensional mapping from the difference between two 
photometric bands with S/N=20. (e) Two-dimensional mapping from the difference between two photometric bands with S/N=100. Plots 
of integrated sensitivity Sj as a function of latitude are attached for respective cases. Note that Sj does not depend on longitude because 
spin rate is much faster than orbital motion. 
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and snow cover fraction (panel (c)). The intensive region 
in the south hemisphere is not seen in the resultant map. 
This is a natural consequence because the region near 
the south pole is insensitive to the data for this specific 
geometry {( = 23°.4, Beq = 90°, i = 45°). The integrated 
sensitivity^ Sj^ to data is plotted as a function of lati- 
tude at the right of each map. The interpretation of the 
reconstructed albedo map should be accompanied with 
the sensitivity plot. 

Surface inhomogeneity is expected to be uncovered by 
mapping the difference between two photometric bands 
since the reflection spectra of clouds are roughly con- 
stant against wavelengths from visible to near infrared 
bands (see Figure 3 in Paper I). In this case, we substi- 
tute dA — dB and ttia — for d and m in Equations 
(p!6|) — ([2T]) where subscripts A and B indicate different 
bands. We again employ the criterion a < 90° for usable 
data because the contribution of clouds is not completely 
compensated by taking difference of two bands. The re- 
constructed maps from the difference between two pho- 
tometric bands with S/N=20 (Figure |4|) are displayed in 
panel (d) in Figure[6l The map for NIR— Blue (left panel) 
can be interpreted as a proxy of the land distribution, 
since most of the land components have higher albedo 
than ocean at the near-infrared range (Figure 3 of Paper 
I). Indeed one can see the bright regions corresponding to 
the Africa continent and other continents (South Amer- 
ica, North America, and Eurasia). On the other hand, 
the map from NIR— Orange (middle panel) is sensitive 
to the vegetated area as well because of the vegetation's 
red edge. The bright regions in the NIR— Orange trace 
the vegetated area of Africa, South America, and Eurasia 
well. 

We emphasize the fact that the SOT maps the albedo 
contrast over planetary surface and the map retrieved by 
the SOT is independent of any assumptions of surface 
components. The interpretation of the albedo contrast is 
in our hands and depends on our knowledge about reflec- 
tion spectra of materials. Hence the degeneracy among 
components with similar reflection spectra is inevitable. 
For instance, distinguishing clouds from snow is diflicult 
from the bands we considered because they both can have 
high and flat reflection spectra in the optical range. 

Panel (e) displays the same results as panel (d) except 
that imposed S/N per data point is 100. The maps with 
S/N= 100 data have higher spatial resolution than those 
with S/N= 20 data and the continental distribution are 
clearly recovered. 

In short summary, the SOT can also be applied to 
Earth-twins in inclined orbits with observations at the 
phases larger than 90°. These phases are also favorable 
from the view point of observation because the intensity 
of planet is larger at smaller phase angle. 

4.3. Measurements of Obliquity and Equinox 

In this section, we test the obliquity measurement in- 
troduced in Section 12.41 Unlike Section 14. 2| we regard 
not only rnj but also ( and 6eq as fltting parameters (but 
the orbital inclination i remains to be flxed). We create 
mock light curves of Earth- twins with varying i, Oeq 
and perform the SOT on them respectively. The same 
data sets for cloud cover, surface reflectivity, and snow 
cover are adopted as input in the simulations though the 
planetary climates may be admittedly inconsistent with 



obliquity when C 7^ 23°. 4. 

Figure [71 shows the measurement of ( and 6eq from 
the NIR— Blue mapping as a function of input value of 
As can be seen, the estimated values are well corre- 
lated to the input values. For the highly oblique planets 
(C ~ 90°), the estimation works remarkably well, con- 
sistent with the face-on cases discussed in Paper I. De- 
pending on the geometry of observation, there are sys- 
tematic biases in some cases. Identifying the exact cause 
of the systematics is complex and we postpone the prob- 
lem of precise measurement of obliquity until a future 
paper. Nevertheless, the fact that the SOT can roughly 
measure the obliquity and equinox is remarkable since 
the obliquity measurement is signiflcantly difficult with 
other method. 

We also tried the same estimation with single-band 
light mapping, but the result is much less reliable. Not 
only latitudinal but also longitudinal inhomogeneity just 
like the land distribution of the Earth appears to be es- 
sential for measuring the planetary spin axis. 

4.4. Prospects with Future Instruments 

Since the development for direct imaging observation 
is currently ongoing and not yet put into practice, we es- 
timate the required aperture size of space telescopes for 
the SOT simply as follows. We consider readout noise, 
dark noise, and noise from exozodiacal light in addition 
to the photon noise of planetary light itself, and calcu- 
late S/N as a function of aperture size with Equations 
(21) -(24) of lKaw^ahara fc Fuiiil pOlOf ). The leakage of 
the light from the host star might be a dominant noise 
source for some instruments but depends on the situa- 
tion. Considering that the occulter concept can safely 
neglect the leakage of the host star, we do not include it 
in our S/N estimation. The input values for the parame- 
ters to specify the instrumental property are summarized 
in Table El 

Based on these assumptions, S/N of 20 with an Earth- 
sized planet at the distance of 10 pc will require ~5 m 
aperture, and that of 100 will be equivalent to ^ 15 m 
aperture size (Figure [8]). This is signiflcantly smaller 
than the scale proposed for future missions aiming to 
obtain spatially-resolved image of exoplanets in a direct 
manner (e.g. hypertelescope; Labeyrie 1999). This is the 
outcome of the SOT which does not resolve the image of 
the planet instantaneously but instead solves the inverse 
problem of annual light curves. 

5. CASE STUDIES 

In this section, we apply the SOT to different conflgu- 
r at ions. We consider different geometries with different 
combinations of orbital inclination and obliquity (Sec- 
tion [5ll]), and different continental distribution based on 
ancient Earth maps (Section 15. 2p . 

5.1. Different Geometry 
5.1.1. Upright Earth- Twin 

For the face-on case examined in Paper I, the two- 
dimensional image can be retrieved for non-zero obliq- 
uity. For inclined {i ^ 0) cases, however, even non- 
oblique planets (i.e., C = 0) can provide two-dimensional 
information of the planetary surface in principle because 
the visible region on the surface changes according to the 
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Fig. 7. — Estimated value of obliquity, C(est) (top row), and simultaneously estimated value of equinox phase, Beq(est) (bottom row), 
as a function of input value of obliquity, C(in)- Different colors show the result with varying input value of equinox phase, Beq. Three 
inclinations, i = 30° (left column), i = 45° (center column), and i = 60° (right column), are studied. Signal-to- noise ratio is 20. 



TABLE 2 

Observation Parameters for Noise Estimation 



Symbol 


Quantity 


Value 


Unit 


^exp 


Exposure time 


24/30x3600 


s 


n 


Folded days 


6 


days 




Sharpness 


0.08 




h 


End-to-end efficiency 


0.5 




V 


Dark rate 


0.001 


counts s"-*^ 


K 


Read noise 


2 


Vcounts read"-*^ 


QE 


Quantum efficiency 


0.91 






Zodiacal light in magnitude 


23 


mag arcsec"^ 


Z 


Zodiacal light 


2 





orbital revolution as shown in Figure [2fc). To demon- 
strate such cases, we simulate light curves of a mock 
Earth with zero obliquity with orbital inclination i = 60° 
(Fig. [2Ic)). In this case, we may ignore the seasonal vari- 
ation of polar snow covers and surface BRDFs (we adopt 
the BRDFs of March all year around) but the same input 
data of clouds are used. 

Figure [9] displays the resultant maps from NIK— Blue 
and NIR— Orange for i = 60° and 0, 



eq 



90°. They 

shows as good correspondence to the real surface inho- 
mogeneity although the resolution is slightly worse than 
the case discussed in Section IH 

5.1.2. Tidally Locked Earth-Twin 

For future projects of direct imaging of HZ planets, 
those around late-type stars are potentially favorable 



targets since late-type stars are most abundant in the 
universe. Since HZs around the late-type stars are so 
close to the star, pla nets there are likely tidally locked 
(|Kasting et al.|[l993f ). While the illuminated area does 
not change for the tidally- locked planet, the visible area 
does change unless a planet is in a face-on orbit (Figure 
[2fc)), which makes the SOT applicable to such cases as 
well, though the reconstructable area is inevitably lim- 
ited to the illuminated half of the surface (see Appendix 
B). 

Here we consider two tidally locked cases (i.e., C = 
0°,cJspin = ^orb)- ouc is an atmosphere- less Earth to see 
the results for the idealized case. Another is an Earth- 
twin with cloud cover same as the previous section. In 
both cases, again, we ignore seasonal variation of snow 
cover and surface albedos. Although those inputs may 
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Fig. 8. — Signal-to-noise ratio as a function of aperture size on 
the assumption of parameters listed in Table [2] 
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Fig. 9. — Two-dimensional mapping from NIR— Blue (left) 
and NIR— Orange (right) light curves in the case of {(", Beq,i} = 
{0°, 90°, 60°}. Data at o; < 90° are used in the analysis. 
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Fig. 10. — Left: albedo mapping from NIR band light curve of 
a hypothetical tidally locked atmosphere-less Earth. Number of 
sampling per orbit (with constant intervals) are 200 and all data 
are used in the analysis (A^data = 200). S/N= 100 is imposed. 
The white dotted lines are boundaries of illuminated and visible 
area. Right: albedo mapping from NIR— Blue light curve of a hy- 
pothetical tidally-locked planet with Earth-twin surface and cloud 
cover. Number of sampling per orbit (with constant intervals) are 
400 and data at o; < 90° are used in the analysis (A^data = 200). 
S/N= 20 is imposed. Assuming S/N= 100 does not improve the 
maps significantly. 

be unrealistic for actual tidally locked Earth analogs, we 



avoid getting deeply involved since the aim of this section 
is to see the basic applicability of the SOT to the tidally 
locked situation. 

We assume that the sampling frequency is 200/orbit 
and 400/orbit, respectively, taking into account the fact 
that the orbit of tidally locked planets is much shorter 
than Earth's (for instance, the orbital period of an Earth- 
mass planet around a star with M = O.54M0 is 2400 
hr). Accordingly, we reduce the number of pixels to 
Mpix = 12-4-4 = 192 (A^side = 4) in order to avoid 
under deter mined problems. 

The left column of Figure [10] displays the recovered 
maps from the atmosphere-less light curve at the NIR 
band with S/N= 100 in the case of z = 60°. For 
demonstrative purpose, we show the maps with ^offset = 
0°, 90°, 180°, and 270°. This inversion is in a sense 
straightforward since our simulation assumes that the 
surface is Lambertian. As expected, the continental dis- 
tribution is recovered fairly well. 

The right column shows the results of mapping from 
the NIR— Blue light curves of a tidally locked Earth- twin 
with cloud cover at a < 90°. Though the resolution 
is now highly deteriorated, there is still inhomogeneity 
which certainly traces from the real features. In this 
case, increasing S/N to 100 does not improve the image 
quality significantly. 

5.2. Different Land Distribution 

In this section, we consider the sensitivity of recon- 
struction to the continental distribution. One of the 
ways to do this is t o use the continental d istribution of 
ancient Earth (see ISanromd &: Pallel[20T2h . We adopt 
the Earth map in the latest Jurassic period (150 Ma) 
and that in the early Camb rian period (540 Ma) taken 
from the paleomap project^ (Scotese"2001) shown by the 
left column of Figure [TTJ We simply assume that every- 
where on the continents has the same soil-like spectrum 
with refiectivity 0.1 in Blue band and 0.3 in NIR band. 
Since the primary features of cloud pattern are controlled 
by the global circulation which depends largely on the 
spin rotation rate and planetary radius, one may expect 
that global cloud pattern is not completely different in 
these periods. Therefore we simply overlaid clouds of the 
present Earth on the ancient continental configuration. 

The middle column of Figure [11] displays two- 
dimensional maps retrieved by the same procedure from 
the NIR-Blue light curves with S/N= 20. The resultant 
map of the Jurassic Earth shows the primary continen- 
tal mass concentrated on the center and a smaller feature 
at the northeast, which trace the input map with lower 
resolution. Likewise, the map of Cambrian Earth also 
recovers the input feature. The difference among the 
resultant maps of Earths at different epochs (Figures [6] 
and [TT]) are evident. These results imply that the SOT 
is applicable to various continental configurations. 

6. SUMMARY AND DISCUSSION 

In this paper, we extended the SOT, a method to re- 
cover two-dimensional information of exoplanetary sur- 
face from annual light curves proposed in Paper I, to be 
applicable to various geometric cases. In Paper I, we 
considered a planet in a face-on orbit, for which global 

^ http://www.scotese.com/ 
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Fig. 11. — Left: input data of land distribution based on paleomap of the latest Jurassic period. Right: reconstructed map from 
NIR— Blue light curves with signal-to-noise ratio 20 based on the latest Jurassic (upper panels) and early Cambrian (lower panels) land 
distribution. Geometric parameters are {C, Beq, i} = {23°. 4, 90°, 45°}. Data at o; < 90° are used in the analysis. The integrated sensitivity 
Sj is plotted for each case in the right column. 



maps with highly oblique planets are available. In this 
paper, we extend the applicable range of this method 
to Earth- twins in an arbitrary inclined orbit, which, in 
some geometric cases, enables the sketch of almost global 
surface map of low-obliquity planet like the Earth. Gen- 
erally, the effect of anisotropic scattering by cloud can 
cause bias, but we have shown that omitting the data at 
crescent phases allows us to reasonably recover the sur- 
face inhomogeneity of Earth- twins. In the same way as 
Paper I, an albedo mapping of an Earth-twin from single- 
band photometry primarily contains the information of 
the distribution of cloud and snow. Information of bare 
surface under the clouds can be estimated by mapping 
the difference of light curves of two different photometric 
bands, thanks to the fact that cloud reflection spectra 
are roughly independent of wavelengths. On the simple 
assumptions for observational noise described in Section 
14. 4[ a space telescope with aperture size ~ 15 m will al- 
low as to obtain two-dimensional map of an Earth-twin 
at a distance of 10 pc with satisfactory spatial resolution 
from continuous 1/2 year observation. 

We tested the reliability of the simultaneous obliq- 
uity measurement with mock light curves of Earth-twins 
with varying directions of the spin axis and the observer. 
We confirmed that the global inhomogeneity of surface 
albedo allows us to reasonably estimate the planetary 
obliquity. 

In the sense that primary features in scattered light 
of the Earth come from clouds because of their global 
cover (typically more than a half of the entire sur- 
face is covered by clouds with optical thickness larger 
than 5), clouds can be regarded as a substantial noise 
source if w e aim t o know the surf ace components 
(iFord et al.ll2 001: Pall e et al.ll2 008: Oa klev &: Cashll2009l : 
ICowan et aO2009. : ..Fuiii et al.lboiO. 201^. It is even an- 



noying since cloud pattern is variable in the timescale of 
^ hour and there is also a seasonal variation. Indeed, the 
short-time variation of cloud cover causes ~10% fluctua- 
tion on diurnal light curves, which is comparable to the 
observational noise we typically imposed (~ 12%). Nev- 
ertheless, it is known that the annual mean of cloud cover 
has a characteristic pattern as shown in the bottom panel 
of Figure [6l which is related to the global atmospheric 
circulation of the Earth. Reconstructed albedo maps are 
likely to infer such averaged features of cloud cover since 
the SOT utilizes long-term light curves. Therefore we are 
no longer bothered by temporarily localized cloud cover. 

The SOT has a potential to reveal the presence and 
rough distributions of the continents. For life on the 
Earth, the presence of continents has played an essen- 
tial role in evolution of life as a source of nutrients by 
weathering. Nutrients such as phosphorus are indispens- 
able to microbes in ocean such as cyanobacteria^ who has 
been creating oxygen in Earth's atmosphere. Therefore, 
discovery of continents on HZ planets may encourage us 
to search for life there. In addition, planets with con- 
tinents can be advantageous targets for the detection of 
direct indicators of photosynthesis, because the signature 
of chlorophyll of microbes in ocean is more difficult to de- 
tect than the red edge of land plants (|Knackell2003l ) . The 
distribution of continents, if exist at all, and the plane- 
tary obliq uity can gs:reatly affect the clime of th e planet 
(e.g., Willi ams K asting 1997: Abe et al.ll2005f ). Know- 
ing these properties via the SOT will potentially enable 
the synthetic discussion of habitability from a variety of 
aspects. 
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APPENDIX 
A. WEIGHT PUNCTION 
Here we present the exphcit form of the weight function in Equation ([9]), 

W{c^, 0- 0) = Wi{^, e- ^; O) Wv{^, 0; ^) 

Wi ((/), 6>; <l>; 6) = max{sin (6 - 6eq) cos ( sin sin + sin (6 - 6eq) sin ( cos 
+ sin 6> cos (6 - 6eq) cos + <l>) , 0} 
Wv ((/), ^; = max{ — sin i sin Oeq cos ( sin sin {(j) -\- ^) — sin i sin Oeq sin ( cos 

— cosi sinC sin 6 sin (0 + <!>) + cosicos^ cos^ + s'mO s'mi cos + <l>) cos6eq,0}. 



The boundary curves of the visible and ihuminated area are expressed as 

^ay((^ + <l>) ± 



<^y,±(^) = 2tan"^ 



<!>) 



(Al) 
(A2) 
(A3) 

(A4) 



ay ((/>) = — COS ( sin i sin (j) sin Oeq — sin ( cos i sin cj) + sin i cos <p cos 0eq 
by {(j)) = (— cos ( sin i sin (j) sin Oeq — sin ( cos i sin + sin i cos cos 0e 

+ (cos ( cos i — sin ( sin i sin 0eq)^ 
cy {(j)) = COS C COS i — sin C sin i sin 6eq, 
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and 



[ COS ^ 

ai {(/)) = 



a/(0 + <!>)] 

cos C sin (9 — 6eq) sin 



for < e - Oeq < TT 
for TT < 6 - Oeq < 27r 



(A5) 



COS (6 — 9eq) COS( 



/sin (6 — 6eq) sin C + (cos(6 — 9eq) cos ^ + cos sin(6 — 6eq) sin ^) 
Using the spherical coordinate {(po^Oo) to the observer on the planetary surface, one can rewrite Equation ()A4p to 

cos {(j) — (j)o) sin^o 



I cos^ + cos^ [(j) — (j)o) sin"^ ^ 

We note that = 27r — ^ and 6>y,±(^) does not depend on 6. 

Converting bq to the spherical surface at ^ = (i.e. Rx{—C)eo) , we obtain the transformation formula to {(po.Oo) 

as 

cos = cos i cos C — sin i sin (" sin 0eq (A7) 
sin sin ((/)o + = — cos i sin — sin i cos ^ sin Oeq (^8) 
sin^o cos ((/)o + ^) = sinz cosGeq. (^9) 

B. RECONSTRUCTABLE AREA 
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Fig. 12. — Contour plot of the cover fraction of the illuminated and visible surface over a planetary year as a function of i and Beq- 
contours show /iv = 0.6,0.7,0.8,0.9, and 0.99. Each panel displays different obliquity corresponding to C = 90°, 60°, 45°, and 30°. 



Solid 



Figure [12] shows the contour of the area cover fraction of the surface sampled at least once during a planetary year 
/iv, that is the maximum reconstructable area, for various geometric configurations. We also derive the cumulative 
probability distribution function p{> fiy). The whole region of the visible area is covered in a planetary year, we 
only have to consider the cover fraction of the visible area {p{> /iv) does not depend on the obliquity). Considering 
diurnal rotation, we can express the cover fraction as a function of Oq^ 

/IV = (Bl) 



Using the uniform distribution of the line of sight p{fiv) df] 



IV 



p{> fiv)= / P{fiv)dfiv = / sini9o 

/l V /l V 



sin6>o dOo and Equation (|B1[) . we obtain 
dOo 



dh 



IV 



dh 



IV 



2(2/iv - 1) 



z#iv = 2v/(l-/iv)/iv. 



(B2) 



//IV Vl-(2/iv-l)^ 

The p{> /iv) does not depend on the obliquity ^ and one can get above 95% (80%) of the surface information in 43.5% 
(80%) for the randomly selected exoplanets, in principle. 
When a planet is in a synchronous orbit, a half of the planetary surface is always dark. We obtain 

1 + sinz 
/IV • 



The CPDF of the tidally locked planet is expressed as 

1 

di 



Ptl(> /iv)= / smz 



dh 



IV 



dhv 

^ v/2(-l + 4/iv) 
hiw \//iv(l - 2/iv) 
Both p{> /iv) & Ptl(> /iv) are shown in Figure [T3l 



dhv = 2^2(1 -2/iv)/iv. 



(B4) 
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Fujii and Kawahara 



P ( > /iv) 




Fig. 13. — Cumulative probability distribution function (CPDF) of the cover fraction /iv on the assumption of the random line of sight 
r2) oc sini), p{> /iv) = 2^ /iv(l — /iv)- This function is independent of the obliquity C- The probability p(/iv > 95%) is 43.5% and 
p(/iv > 80%) is 80%. The dashed line indicates the CPDF for tidally locked planets: ptl(> /iv) = 2v^2/iv(l - 2/iv). The probability 
Ptl(/iv > 47.5%) is 43.5% and J9tl(/iv > 40%) is 80%. 



C. BAYESIAN INTERPRETATION OF THE TIKHONOV REGULARIZATION 

Following iTarantolal (|2QQ5l ). we discuss the Bayesian interpretation of the Tikhonov regularization. The probability 
of a posterior information in the model space is given by 

PM{m) oc pM{m)pD{g{m)), (CI) 

Here we neglect modelization uncertainties. Introducing a Gaussian model for the a priori information on the model 
parameters with the covariance matrix Cm, 



pM{m) = exp 

/(2^)^,detCM 



(m — rfi)^Cy^{m — rfi) 



and assuming the statistical noise is an independent Gaussian noise with the covariance matrix Cd 



PD{d) = I exp 

/(2^)^,,,detCD 



one can obtain 



PM{m) (X exp [e'^Cjj^e + (m - rh)^Cyl{m - rh)] |. 



(C2) 



(C3) 



(C4) 



The model which maximizes a posteriori probability is given by minimizing 

Q = s^C^^s + (m - rhfC^^im - rh). (C5) 
If taking the independent Gaussian noise (Cojij = cr|% and (CM)ij = X~'^5ij^ we obtain the Tikhonov regularization 



fc^I^ + A^lm-mp. (C6) 



D. THE L-CURVE CRITERION 
The regularization parameter A is chosen based on "the L-curve criterion" following iHanse ^ (|2QTqI) . The ij-curve IS 

a parametric plot of the model norm |mest,A — ''^1 and the residuals |(5mest,A — d| as a parametric function of A. In 
the left panel of Figure [lH we demonstrate the L-curve of the mock observations of an atmosphere-less Earth with 
C = 90°, Oeq = 0°,i = 0° and 3% noise. As A increases, |mest,A — ''^1 rapidly decreases until B (the L-curve's corner), 
and then the slope becomes gentler. The large value of |mest,A — ''^1 (small |(5mest,A — ^|) implies that the perturbation 
error due to observational noise is large, while large |(5mest,A — d\ (small |mest,A — ''^D, indicates that the bias due 
to the regularization (regularization error) is large. The L-curve corner serves as the point of good balance between 
th ese errors. Th e L-curve corner is defined by the maximum curvature in the log- log space, as described in Appendix 
E (|Hansenll2QlQf ). Adopting different values for A to Equation (|2T]) , we see the A dependence of the recovered models 
mest,A m the right panels. As shown in the left panel, the L-curve criterion gives A = 10^"^^ at the maximum curvature 
point (B). Right panel of Figure [14] shows the recovered maps with A = 10~^"^ (A), the L-curve corner A = 10^"^^ (B), 
and A = 10^-^ (C). The recovered map for A exhibits large perturbation noise, while the resolution of the model is 
degraded for the case (C). The L-curve criterion balances these two opposite effects (B). 
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Fig. 14. — A dependence of the best-fit model. The left panel shows the L-curve for the mock observations with C, = 90°, Beq = 0°, i = 0°, 
and the additional 5% noise. The input data are a land/ocean binary map of the Earth, which has 1 on land and on ocean. A solid curve 
indicates the parametric plot of |Gmest,A ~ ^1 cii^d |mest,A ~ ''^l cis a function of the regularization parameter A. A dotted curve shows 
the curvature of the L-curve, c(A). The definition of c(A) is written in Appendix E. We took three different regularization parameters: 
A = 10~°-^ (A), 10°-^^, which corresponds to the maximum curvature (B), and A = 10-*^ (C). The recovered models are displayed in the 
right panels: the upper, middle, and lower panels correspond to A, B, and C. 

E. MAXIMUM CURVATURE OF THE L-CURVE 
Here, we briefly describe the derivation of the curvature of the L-curve following iHanse 3 (|2Q1Q[ ). We first introduce 
^ = |mest,A — and p = |Gmest,A — The curvature of the L-curve is defined in the log|Gmest,A — epi- 
log |mest,A - rh\ space, 

^ (logp)^(logO^^-(logp)^^(logO^ 

^ ^" [(iogpy2 + (iogO'2]3/2 ' 

where the prime indicates the derivative by A. The derivatives of ^ = |mest,A — iti\'^ and p = |Gmest,A — are 



A 



P = 

Wi{X) = 

Using Equations ()E2|) and (jE2|) . one obtain 

c(A) 



A^datE 

^2 



GjkiTik 



(El) 

(E2) 

(E3) 

(E4) 
(E5) 

(E6) 



